Free expansion of Bose-Einstein condensates with quantized vortices 
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The expansion of Bose-Einstein condensates with quan- 
tized vortices is studied by solving numerically the time- 
dependent Gross-Pitaevskii equation at zero temperature. 
For a condensate initially trapped in a spherical harmonic 
potential, we confirm previous results obtained by means of 
variational methods showing that, after releasing the trap, 
the vortex core expands faster than the radius of the atomic 
cloud. This could make the detection of vortices feasible, by 
observing the depletion of the density along the axis of rota- 
tion. We find that this effect is significantly enhanced in the 
case of anisotropic disc-shaped traps. The results obtained 
as a function of the anisotropy of the initial configuration are 
compared with the analytic solution for a noninteracting gas 
in 3D as well as with the scaling law predicted for an inter- 
acting gas in 2D. 

03.75.Fi, 03.75.-b , 05.30.Jp 



I. INTRODUCTION 

Since the discovery of Bose-Einstein condensation in 
trapped gases of alkali-metal atoms JjJ, one of the pri- 
mary goals of both the experimental and theoretical ac- 
tivity has been the search for superfLuid effects. A clean 
connection among Bose-Einstein condensation and su- 
perfluidity would be, in fact, of fundamental importance 
from a conceptual viewpoint, for it would improve our 
understanding of the properties of quantum many-body 
systems. In this context, attempts have been made to 
produce and detect quantized vortices in these trapped 
gases Q and experiments are currently running in dif- 
ferent laboratories. At the same time, several theoreti- 
cal papers have been written on the expected features of 
these vortices (see, for instance, Refs. and references 
therein). 

One of the difficulties in observing a vortex in a 
trapped condensate is that its core, i.e., the region where 
the density is depleted by the centrifugal force associated 
with the quantized vorticity, has a radius of the order of 
the healing length £. This quantity is significantly smaller 
than the radius of the cloud, especially when the number 
of atoms, N, is large. This makes a direct observation 
rather difficult, due to the finite resolution of the optical 
imaging devices. 

Among the various methods suggested for detecting 
vortices, one consists in letting the condensate expand 



freely by switching-off the trapping potential. It has re- 
cently been shown |?J that the core of the vortex may 
expand faster than the radius of the cloud, making the 
observation of the corresponding density depletion feasi- 
ble. In order to give quantitative estimates of this effect, 
in Ref. Jjj the Gross-Pitaevskii equation has been solved 
by using an approximate trial wave function and the vor- 
tex core size has been studied in the case of a gas initially 
confined in a spherical trap. 

In the present work we study the same problem, 
namely the dynamics of the expanding condensate with 
a vortex inside, by exactly solving the Gross-Pitaevskii 
equation and exploring also the effects of the anisotropy 
of the initial configurations. Our calculations for the 
spherical case closely agree with the results of Ref. 0, 
showing that the magnification of the core size is not an 
artifact of the variational wave function used by those au- 
thors. In the case of anisotropic configurations, we find 
different behaviors depending on the form of the trap. 
When the initial confinement is weaker along the axis 
of rotation (cigar-shaped traps) the axial expansion is 
slower, while the radial one is similar to that of a two- 
dimensional condensate: the core size and the radius of 
the cloud increase with almost the same speed, follow- 
ing the scaling behavior expected for the solutions of the 
Gross-Pitaevskii equation in 2D Conversely, when 
the initial confinement is stronger along the axis of ro- 
tation (disc-shaped traps), the axial expansion is faster 
and the core size is found to increase much more than the 
radius of the cloud in the first instants of motion. The 
ratio between the radii of the core and the cloud starts 
from a rather small value, which depends on N and is 
close to the prediction for a 2D condensate, and then in- 
creases quickly to the value 0.282, which is the analytic 
result for an ideal gas in 3D. All these results will be dis- 
cussed in details in section III, while in the next section 
we will introduce the basic formalism and the numerical 
procedure. 



II. SOLVING THE GROSS-PITAEVSKII 
EQUATION 

The dynamics of a dilute and weakly interacting Bose- 
Einsten condensate at zero temperature is well described 
by a mean-field theory, in which the motion of all the 
atoms is assumed to be determined by a single complex 
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function <I>(r,t), also called order parameter, playing the 
role of a macroscopic wave function and obeying the fol- 
lowing equation: 



h 2 v 2 

2m 



*(r,t) 



(1) 



This is known as Gross-Pitaevskii (GP) equation ||. The 
mean-field interaction enters through the term propor- 
tional to the particle density n(v, t) — |$(r, t)\ 2 , the cou- 
pling constant g being related to the s-wave scattering 
length a by g = Air% 2 a/m. For an axially symmetric trap 
the confining potential V ex t (r) can be written in the form 
V cxt (p, z) = {m/2)[uj 2 p p 2 + uj 2 z z 2 ], with p = {x 2 + y 2 ) 1 / 2 . 
The ratio between the axial (lo z ) and radial (lo p ) fre- 
quencies, A = uj z /ujp, is a useful parameter fixing the 
anisotropy of the trap. The derivation and several ap- 
plications of the GP equation (fTt) are reviewed in Ref. 
1- 

In general the order parameter entering the GP equa- 
tion (Q) can be written in the form 



exp[iS(r, t)] 



(2) 



where the phase S(r, t) can be used to define a velocity 
field through 



v(r,t) = -VS(r,t) 



(3) 



When a quantized vortex is present in the condensate, 
with its axis of rotation along z, the atoms flow with 
tangential velocity v = nh/(mp), each one having an- 
gular momentum L z = nh, where k is the quantum of 
circulation. This means that, from Eq. (j^), the order pa- 
rameter has to depend on the angle 0, around the z-axis, 
in a simple way: 



3>(r, t) = z, t) cxp[iK(, 



(4) 



Inserting this expression back into the GP equation ([!]), 
one gets 



h 2 v 2 

2m 
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2m p 



ih^(p,z,t) = 

+j(uy + uj 2 z 2 )+g\<f(p,z,t)\ 2 ] 9(p,z,t), (5) 

which corresponds to the GP equation with a centrifugal 
force included. Due this centrifugal term, the solution 
ty(p,z,t) must vanish on the z-axis. The region where 
the density is depleted, close to the vortex line, is named 
vortex core. In a stationary state, its radius is of the order 
of the healing length £ fixed by the balance between the 
mean-field energy and the kinetic energy associated with 
the gradient of the density (quantum pressure). In a uni- 
form gas, the healing length is given by £ = (87rna) -1 / 2 , 
where n is the density. One can use the same expression 



in order to estimate the core size in trapped gases by 
taking for n the central density of the condensate with- 
out vortices. With the parameters of the current exper- 
iments the predicted core radius turns out to be very 
small indeed, making the direct observation of a vortex 
rather difficult with the presently available optical de- 
vices. However, when the condensate is let to expand 
by switching-off the confining potential, the vortex core 
expands as well. We want to describe this process by 
numerically solving Eq. (|5|). 

The GP equation (|5|) can be formally rewritten in this 
way: 



ih- 



dt 



(6) 



where the effective Hamiltonian H is the operator en- 
closed in square brackets in (^). Equation (^) has the 
form of a nonlinear Schrodinger equation, since H con- 
tains the unknown \&. There exist several techniques to 
solve it numerically. One of them consists in propagating 
the wave function in time using small time steps At in 
order to approximate the time evolution of \& with the 
equation 



iAt \ ( iAt 



(7) 



where H and ^> n are the effective Hamiltonian and the 
order parameter, respectively, both calculated at the n- 
iteration. Actually, since the diffusion of ^> occurs in 
both the axial and radial directions, one splits each time 
interval in two parts, one for the evolution along z and 
the other along p. In order to propagate \& in each di- 
rection, we separate the effective Hamiltonian into axial 
and radial components, equally dividing the mean-field 
potential (this separation is somewhat arbitrary and does 
not affect the final results). Equation for both the 
axial and radial motion at each time step, can be solved 
by mapping the order parameter on a two-dimensional 
grid of points, N p x N z , in such a way that its solution 
becomes equivalent to a matrix diagonalization. This is 
known as Crank-Nicholson differencin g m ethod with al- 
ternating direction implicit algorithm |10[ and was used 
for the expansion of a trapped condensate first by Hol- 
land and Cooper JO]. 

For a vortex in a trap with radial frequency uj p , the rel- 
evant length scale is a p — [h/^miOp)] 1 ' 2 . Available con- 
densates have radii of the order of 1 to 10 a p , depending 
on the strength of the mean-field interaction fixed by the 
parameter Na/a p 0). We consider condensates with re- 
pulsive interaction (a > 0) and with Na/a p from to 
20. We choose a grid having typically 100 x 200 points 
with spacing smaller than 0.2a p , which is sufficient for 
an accurate description of the density profile of the con- 
densate including the vortex core. The derivatives in the 
effective Hamiltonian H are expressed through finite dif- 
ference 5-points formulas p2j. 
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The initial configuration at t = corresponds to the 
stationary order parameter of the condensate in the trap 
for a given k. It can be obtained by using the same al- 
gorithm, but diffusing the function W in imaginary time 
instead of real time (t — > — it). One starts from a rea- 
sonable trial wave function and let W converge to the 
stationary state. The same can be done by using an ex- 
plicit differencing algorithm instead of the implicit one 
based on Eq. (]?]). The explicit method has been used, 
for instance, in Ref. fl3|| and turns out to be much faster. 
We used both methods as a test of the numerical code. 

For t > the confining potential in Eq. (||) is removed 
and hence the function \l/ expands freely, subject only to 
the mean-field repulsion and quantum pressure. A CPU 
time of the order of 2 hours on a workstation is typically 
needed to follow the expansion of the condensate for a 
time of the order of lOw" 1 , which is enough to see the 
effects that we are going to discuss. 

The expansion of an ideal gas (a = 0) is used as one of 
the tests we made on the numerics. In this case, one has 
simple analytic results. The t = configurations are the 
eigenstates of the harmonic oscillator in axial symmetry. 
The k = case is a Gaussian, while the vortex state 
k = 1 corresponds to the eigenfunction with azimuthal 
angular momentum m = 1. When the trap is released, 
their expansion is just the dispersion of free wave packets. 
We checked that our numerical integration gives results 
practically indistinguishable from the analytic ones. In 
order to verify the proper inclusion of the mean-field in- 
teraction in the code, we made the calculation for the ex- 
pansion of a condensate without vortices, using the same 
parameters of a previous work [fl4[ and finding exactly 
the same results. As a further check, we looked at the 
behavior of the energy of the condensate as a function of 
time. It is well known that the Crank-Nicholson method 
may suffer from a possible amplification of random nu- 
merical noise, due to nonlinearity, leading to a spurious 
input of kinetic energy and limiting the stability of the 
algorithm to short time. However, we checked that this 
time is long enough to ensure the stability of the solutions 
in the time intervals here considered. 



III. RESULTS 

First we give two examples of expanding vortices. In 
Fig. |l| and || we show a sequence of contour plots of 
the density n(p,z,t) — \^(p,z,t)\ 2 in three instants of 
time (t = 0, 1.5, Slo^ 1 ) and for two different condensates 
(A = 0.2 and 2.5, respectively). The darkest color repre- 
sents the region where the density is larger than 1/2 of 
the peak density at t = 0, and each contour line is drawn 
where the density falls by a factor 2. Both condensates 
have Na/cip — 20, so that the effects of the mean-field 
interaction on the vortex structure are of the same order. 
The empty region close to p — is the vortex core, whose 
radius grows when the atomic cloud expands. The con- 



densate in Fig. is initially confined in a cigar-shaped 
trap. When the trap is released, its motion is mainly in 
the radial direction, the axial size increasing much more 
slowly Jl5|]. The opposite happens for the condensate in 
Fig. 0: the initial configuration is oblate (disc-shaped) 
and the axial motion is faster. 

The shape of the condensate is better seen by plotting 
the column density, that is the density integrated along 
the z-direction, n co \{p,t) = J dzn(p, z,t). This quantity 
is directly measured in the experiments. In Fig. ^, we plot 
n co i for the same expanding condensate shown in Fig. |^. 
To give an idea of the actual size, let us suppose to have 
a trap with frequency lu p = 27r(20Hz) and condensates of 
87 Rb or Na atoms; then the length scale a p is 2.4 pm and 
4.7 /im, respectively. In the same cases, the parameters 
Na/a p — 20 and A = 2.5, as in Fig. |3[ would correspond 
to N ~ 8000 for 87 Rb and N ~ 34000 for Na. 

By using the column density one can give suitable def- 
initions of the radii of both the atomic cloud and the 
vortex core. The former can be identified with the root 
mean square radius, p rms , obtained by integrating p 2 over 
the column density distribution, while the core radius, p Cl 
can be chosen as the value of p where the column density 
first reaches e _1 times the peak value. It is then useful to 
study the evolution in time of the ratio p c j p rms , in order 
to have an estimate of the size of the "hole" made by a 
vortex in typical condensates of current experiments, as 
done in Ref. @ . In Fig. || we show this ratio for different 
values of A and Na/a p . This figure contains the main 
results of the present work. 

First, the horizontal dot-dashed line is the analytic re- 
sult for the free dispersion of a wave packet of nonin- 
tcracting particles. In this case, the motion in the two 
directions is decoupled and the ratio p c /prms is approx- 
imately 0.282 independently of A. The points close to 
this value correspond to the numerical solution of the GP 
equation (|B|) without the mean- field term (g — 0). The 
fluctuations at short time are not due to random noise 
in the evolution of 'J', which remains indeed very close to 
the analytic solution at any time, but to the smallness of 
the vortex core at the beginning of the expansion: both 
the position of the peak and that of p c are located in be- 
tween grid points and the use of interpolation formulas 
yields small errors. As one can see from the figure, these 
fluctuations are rather small and do not affect the main 
results of this work. 

The results obtained by expanding condensates with 
Na/a p = 5 and 20 are also shown. At t = 0, each se- 
quence of points starts from the value of p c /p r ms pre- 
dicted for a stationary state in the trap. This value is a 
decreasing function of both the parameters Na/a p and A. 
An estimate of this dependence can be obtained in the 
so-called Thomas- Fermi limit Na/a p 3> 1, when quan- 
tum pressure is negligible compared to the mean-field 
repulsion. In this case, one finds p c J prms ~ £//°rms cx 
{lh\Na/a p )~ 2 l b . Thus the portion of condensate where 
the density is depleted by the vortex becomes smaller 
and smaller when N and/or A increase. This is why in 
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current experiments the vortex core might be hardly ob- 
servable by simply looking at the density distribution in 
the trap. 

At t > the ratio p c /p r ms increases in different ways 
depending on the anisotropy of the trap. In the case of 
a spherical trap (A = 1), we compare the numerical solu- 
tion of the GP equation (points) with the results of the 
variational calculation by Lundh, Pethick and Smith |?]] 
(solid lines). The agreement is remarkable. The ratio 
Pel Prms is shown to increase rather quickly towards the 
analytic result for the noninteracting gas. As discussed 
in Ref. 0, this is due to the fact that, in the first instants 
of motion, the core size adjust almost istantaneously to 
the value £ = {Sima^ 1 / 2 , with n given by the local den- 
sity just outside the vortex core. This density decreases 
rapidly, yielding an increase of p c faster than the one of 
p rms . After a characteristic decoupling time, this mean- 
field effect tends to vanish, since the gas becomes more 
and more dilute, and the condensate expands as an ideal 
gas, with the axial and radial motions decoupled. The 
ratio Pc/ Prms thus remains almost frozen to its value at 
the decoupling time. 

As one can see in Fig. ||, the magnification of the core 
size depends on both Na/a p and A. In order to em- 
phasize the dependence on the anisotropy of the trap, in 
Fig. H the results for Na/a p = 20 and different values 
of A are plotted together. The lowest dashed line is the 
result we obtain by solving the GP equation in a purely 
2D case with the same Na/a p . In this case, one can 
prove that the dynamics is governed by an exact scaling 
law ||, that is, the condensate preserves its shape with 
a time dilatation of length scales; thus p c / p Tm s remains 
constant. Our results show that the ratio p c j 1 p vms is al- 
most constant even for condensates initially confined in 
prolate (cigar-shaped) traps, the value at t = being 
fixed by the local density near the vortex in the station- 
ary state. This can be understood by noticing that, for 
A « 1, the axial motion is so slow that it does not af- 
fect significantly the behavior of the radial expansion. 
In other words, the lowering of the local density, which 
causes p c to increase, is almost completely determined by 
the radial expansion; the condensate can be thought as 
made by thin slices, orthogonal to the z axis, each one 
expanding as a 2D system. The opposite happens in the 
case of disc-shaped traps (A » 1), where the initial value 
of p c / p rms is indeed close to the one of a 2D system, but 
quickly approaches the 3D ideal gas prediction. This is 
again a consequence of the axial motion. In this case in 
fact, the rapid decrease of the local density around the 
vortex is almost entirely due to the fast axial expansion 
and hence the core size increases much faster than in a 
2D system. 

The fact that the 2D scaling behavior applies to the 
expansion of cigar-shaped condensates rather than the 
disc-shaped ones was not obvious so far. For instance, in 
a recent paper by Castin and Dum (jlfj), focused on the 
properties of vortices in disc-shaped traps, the expansion 
of the condensate has been treated by keeping the con- 



finement along z constant and using the 2D scaling law 
for the radial motion. However, as shown in the present 
work, the motion along z in the actual 3D experiments is 
not expected to be a small correction to the 2D scaling 
law. On the contrary, it is expected to dominate the dy- 
namics of the expansion, yielding values of /9 c /p rms close 
to 0.282 in short time and almost independently of N. 
This may help the direct detection of quantized vortices 
in these inhomogeneous Bose-Einstein condensates, the 
size of the "hole" in the density distribution becoming 
easily larger than the spatial resolution of the imaging 
devices. In order to show how this expanding holes might 
appear in an experimental observation, the column den- 
sity of Fig. H is plotted again in Fig. ^| as a 3D plot, the 
three frames corresponding to t = 0, 1.5 and 3W" 1 . 
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FIG. 1. Contour plot of the density n(p, z, t) = \^f(p, z, t) 2 
in three instants of time (from left to right, t = 0, 1.5, 3u> p ) 
for a condensate with A = luz/uj p = 0.2 and Na/a p — 20 
and having a quantized vortex (ac = 1) along the z-axis. The 
darkest color represents the region where the density is larger 
than 1/2 of the peak density at t = 0, and each contour line is 
drawn where the density falls by a factor 2. The co-ordinates 
z and p are given in units of a p — ^/(muip)] 1 ' 2 , where w p is 
the radial trapping frequency at t = 0. 



FIG. 2. Same as in Fig. [l] but for a condensate in a trap 
with A = ujz/ojp = 2.5. 



FIG. 3. Column density n(p, t) obtained by integrating 
along z the density plotted in Fig. H. The three curves corre- 
spond to t = 0, 1.5, Scjp 1 . 



FIG. 4. Ratio between the vortex core radius and the 
root mean square radius of the condensate, p c /p rm s, as a 
function of time, for three different initial configurations 
(A = 0.2, 1, 2.5) and for Na/a p = 0, 5, and 20. Time is in 
units of uip 1 . Points are the numerical results of the present 
work. The horizontal dot-dashed line corresponds to the ex- 
pansion of an ideal gas in the state m = 1. The two solid 
lines for Na/a p = 5 and 20 in the spherical case, A = 1, are 
the results of recent variational calculations M. 



FIG. 5. Ratio p c /prms as a function of time, for three 
different initial configurations (A = 0.2, 1, 2.5) but the same 
value of Na/a p = 20. The dot-dashed line is the analytic 
result for an expanding ideal gas. The dashed line is the 
result obtained by solving the GP equation in a purely 2D 
system. 



FIG. 6. The column density of Fig. ^ is shown here as a 
3D plot in order to emphasize the structure of the expanding 
vortex. The three frames correspond again to t = 0, 1.5 and 
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